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Abstract — Key questions that scientists and engineers typically want to address 
can be formulated in terms of predictive science. Questions such as: "How 
well does my computational model represent reality?", "What are the most 
important parameters in the problem?", and "What is the best next experiment 
to perform?" are fundamental in solving scientific problems, mystic is a 
framework for massively-parallel optimization and rigorous sensitivity analysis 
that enables these motivating questions to be addressed quantitatively as 
global optimization problems. Often realistic physics, engineering, and materials 
models may have hundreds of input parameters, hundreds of constraints, and 
may require execution times of seconds or longer. In more extreme cases, 
realistic models may be multi-scale, and require the use of high-performance 
computing clusters for their evaluation. Predictive calculations, formulated as 
a global optimization over a potential surface in design parameter space, may 
require an already prohibitively large simulation to be performed hundreds, if 
not thousands, of times. The need to prepare, schedule, and monitor thousands 
of model evaluations, and dynamically explore and analyze results, is a chal- 
lenging problem that requires a software infrastructure capable of distributing 
and managing computations on large-scale heterogeneous resources. In this 
paper, we present the design behind an optimization framework, and also 
a framework for heterogeneous computing, that when utilized together, can 
make computationally intractable sensitivity and optimization problems much 
more tractable. The optimization framework provides global search algorithms 
that have been extended to parallel, where evaluations of the model can be 
distributed to appropriate large-scale resources, while the optimizer centrally 
manages their interactions and navigates the objective function. New methods 
have been developed for imposing and solving constraints that aid in reducing 
the size and complexity of the optimization problem. Additionally, new algorithms 
have been developed that launch multiple optimizers in parallel, thus allowing 
highly efficient local search algorithms to provide fast global optimization. In this 
way, parallelism in optimization also can allow us to not only find global minima, 
but to simultaneously find all local minima and transition points ~ thus providing 
a much more efficient means of mapping out a potential energy surface. 

Index Terms — predictive science, optimization, uncertainty quantification, verifi- 
cation, validation, sensitivity analysis, parallel computing, distributed computing, 
heterogeneous computing 



Introduction 

Recently, a unified mathematical framework for the rigorous 
construction and solution of uncertainty quantification (UQ) 
problems was formulated [OSS 11]. This framework, called 
Optimal Uncertainty Quantification (OUQ), is based on the 
observation that, given a set of assumptions and information 
about the problem, there exist optimal bounds on the uncer- 
tainties. These bounds are obtained as extreme values of well- 
defined optimization problems that correspond to extremizing 
probabilities of failure subject to the constraints imposed by 
scenarios compatible with the information set. 
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An accompanying software framework that implements 
these rigorous UQ/OUQ methods is now posed. 

A rigorous quantification of uncertainty can easily require 
several thousands of model evaluations f{x). For all but the 
smallest of models, this requires significant clock time ~ a 
model requiring 1 minute of clock time evaluated 10,000 times 
in a global optimization will take 10,000 minutes (~ 7 days) 
with a standard optimizer. Furthermore, realistic models are 
often high-dimensional, highly-constrained, and may require 
several hours to days even when run on a parallel computer 
cluster. For studies of this size or larger to be feasible, a 
fundamental shift in how we build optimization algorithms 
is required. The need to provide support for parallel and dis- 
tributed computing at the lowest level ~ within the optimiza- 
tion algorithm ~ is clear. Standard optimization algorithms 
must be extended to parallel. The need for new massively- 
parallel optimization algorithms is also clear. If these parallel 
optimizers are not also seamlessly extensible to distributed and 
heterogeneous computing, then the scope of problems that can 
be addressed will be severely limited. 

While several robust optimization packages exist [JOPOl, 
KROOO], there are very few that provide massively -parallel 
optimization [BMMIO, EKL02, MAT09] ~ the most notable 
effort being DAKOTA [DAKOT], which also includes meth- 
ods for uncertainty quantification [DAKUQ]. A rethinking of 
optimization algorithms, from the ground up, is required to 
dramatically lower the barrier to massively -parallel optimiza- 
tion and rigorous uncertainty quantification. The construction 
and tight integration of a framework for heterogeneous paral- 
lel computing is required to support such optimizations on 
realistic models. The goal should be to enable widespread 
availablility of these tools to scientists and engineers in all 
fields. 

Several of the component pieces of such a framework 
for predictive science already exist, while a few key pieces 
must be constructed ~ furthermore, these packages must 
then be assembled and integrated. Python [GVRPY] is a 
natural integration environment, and is one that readily sup- 
ports the dynamic nature of working across heterogeneous 
resources. By requiring this framework be pure-Python, many 
of the barriers to running on a new platform are removed, 
multiprocessing [MPROC], mpi4py [MPI4P], and pp 
[VVPPP] are selected for communication mechanisms, both 
due to their high level of feature coverage and their relative 
ease of installation. NumPy [NUMPY] is used for algorithmic 
efficiency, and SymPy [SYMll] is used to provide an alternate 
interface for building constraints. Many of the optimization 
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algorithms leverage SciPy [JOPOl]; however like the use 
of Matplotlib [MATPL] for plotting, SciPy is an optional 
dependency. 

This paper will discuss the modifications to the mystic 
[MHA09] optimization framework required to provide a sim- 
ple interface to massively parallel optimization, and also to 
the pathos [MBA 10] framework for staging and launching 
optimizations on heterogeneous resources. These efforts lever- 
age pyre [MAGAl] ~ an component integration framework 
for parallel computing, which has recently been extended to 
distributed communication and management with hydra (part 
of this development effort). This paper will also overview 
a new mathematical framework [OSSll, ALLll, KLLll, 
LOO08] for the quantification of uncertainties, which pro- 
vides a formulation of UQ problems as global optimization 
problems. 

Rigorous Uncertainty Quantification 

Following [LOO08], we specifically take a certification point 
of view of uncertainty quantification. For definiteness, we 
consider systems whose operation can be described in terms 
of N scalar performance measures (7i , . . . , F^v) = 7 G M^. The 
response of the system is taken as stochastic due to the intristic 
randomness of the system, or randomness in the input parame- 
ters defining the operation of the system, or both. Suppose that 
the outcome 7 G A constitutes a satisfactory outcome for the 
system of interest, for some prescribed measureable admissible 
set A C R^. Hence, we are interested in determining the 
probability of failure (PoF) P[7 G A"^]. 

Evidently, for an upper bound to be useful, it must also 
be tight (i.e. it must be close to the actual PoF of the 
system) and accessible by some combination of laboratory and 
computational means. In [ALLll, KLLll], a methodology 
for a rigorous determination of tight upper bounds on the 
probability of failure for complex systems is presented, and 
is summarized below. 

We consider a response function Y = F(X) that maps con- 
trollable system inputs X to performance measures 7, and re- 
lies on a probability of failure (PoF) upper bounds of the con- 
centration of measure (CoM) type [BBL04, LEDOl, MCD89]. 
If McDiarmid's inequality [MCD89] (i.e. the bounded dif- 
ferences inequality) is used to bound PoF, the system may 
then be certified on the sole knowledge of ranges of its input 
parameters ~ without a priori knowledge of their probability 
distributions, its mean performance E[7] = M and a certain 
measure Dq = U of the spread of the response, known 
as system diameter, which provides a rigorous quantitative 
measure of the uncertainty in the response of the system. 

A model is regarded as Y = F{X) that approximates the 
response Y = G{X) of the system. An upper bound on the 
system diameter ~ and thus on the uncertainty in the response 
of the system ~ then follows from the triangle inequality Do < 
Dp -\-Dg-f, and U = Dp -\-Dg-f can be taken as a new - 
- and conservative ~ measure of system uncertainty. In this 
approach, the total uncertainty of the system is the sum of 
the predicted uncertainty (i.e. the variability in performance 
predicted by the model as quantified by the model diameter 



Dp), and the modeling- error uncertainty (i.e. the discrepancy 
between model prediction and experiment as quantified by the 
modeling-error diameter Dq-f- 

In [LOO08], PoF upper bounds of the CoM type were for- 
mulated by recourse to McDiarmid's inequality. In its simplest 
version, this inequality pertains to a system characterized by 
N real random inputs X = (Xi , . . . G £" C and a single 
real performance measure 7 G R. Suppose that the function 
G : R^ R describes the response function of the system. 
Suppose that the system fails when Y < a, where <2 is a 
threshold for the safe operation of the system. Then, a direct 
application of McDiarmid's inequality gives the following 
upper bound on the PoF of the system: 

P[G<a]<expf-2-^j (1) 

where 

M=(E[G]-^z)+ (2) 
is the design margin and 

U = Dg (3) 

is the system uncertainty. In (3), Dq is the diameter of the 
response function. From (1) it follows that the system is 
certified if 




where e is the PoF tolerance, or, equivalently, if 

where CF is the confidence factor. In writing (2) and subse- 
quently, we use the function x+ := max(0,x). We see from 
the preceding expressions that McDiarmid's inequality sup- 
plies rigorous quantitative definitions of design margin and 
system uncertainty. In particular, the latter is measured by 
system diameter Dg, which measures the largest deviation in 
performance resulting from arbitrarily large perturbations of 
one input parameter at a time. Within this simple framework, 
rigorous certification is achieved by the determination of two- 
-and only two—quantities: the mean performance E[G] and the 
system diameter Dq. 

McDiarmid's inequality is a result in probability theory that 
provides an upper bound on the probability that the value 
of a function depending on multiple independent random 
variables deviates from its expected value. A central device 
in McDiarmid's inequality is the diameter of a function. We 
begin by recalling that the oscillation osc{f^E) of a real 
function f : E over a set £" G 7? is 

osc(/,£) = sup{\f{y)-f{x)\ : x,y G E} (5) 

Thus, osc(/,£^) measures the spread of values of / that may be 
obtained by allowing the independent variables to range over 
its entire domain of definition. For functions f : E C R^ R 
of several real values, component- wise suboscillations can be 
defined as 

osCi{f,E) = sup{\f{y) - f{x)\ : x,y G E, Xj=yj for j 7^ /} 

(6) 
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Thus osc; (/,£") measures the maximum oscillation among all 
one-dimensional fibers in the direction of the ith coordinate. 
The diameter D{f^E) of the function / : £" ^ R is obtained as 
the root-mean square of its component- wise suboscillations: 

D{f^E)=\j;^oscj{f,E)j (7) 

and it provides a measure of the spread of the range of 
the function. Thus (6) also us to regard osC;(/,£) as a 
subdiameter of the system corresponding to variable X^, where 
the subdiameter can be regarded as a measure of uncertainty 
contributed by the variable Xi to the total uncertainty of the 
system. 

The attractiveness of the McDiarmid CoM approach to UQ 
relies on the requirement of tractable information on response 
functions (sub-diameters) and measures (independence and 
mean response). Above, it is described how to "plug" this in- 
formation into McDiarmid' s concentration inequality to obtain 
an upper bound on probabilies of deviation. One may wonder 
if it is possible to obtain an "optimal" concentration inequality, 
especially when the available information may not necessarily 
be sub-diameters and mean values. A general mathematical 
framework for optimally quantifying uncertainties based only 
on available information has been proposed [OSS 11], and will 
be summarized here. Assume, for instance, that one wants to 
certify that 

P[G >a]<£ (8) 

based on the information that osc/(G,£') < Z)/, X = 
(Xi, . . . E[G] < and that the inputs Xi are independent 
under P. In this situation, the optimal upper bound ^ {s^md) 
on the PoF P[G > a] is the solution of the following optimiza- 
tion problem 



^(=^md)= sup ^[f{X)>a] 

{f,!i)e^MD 



(9) 



subject to constraints pro vied by the information set 



MD ■ 



Id e J{(Ei)®---®Jl(E^), 
oscK/,^)<A- 



(10) 



where M{E]^ denotes the set of measures of probability on 
E]^. Hence, McDiarmid's inequality is the statement that 

^(^md) < exp (11) 

Similarly, for any other set of information , we have an 
optimal (i.e.) least upper bound on the probability of deviation 



^(j^)= sup \i\f(X)>a\ 



(12) 



The idea is that in practical applications, the available infor- 
mation does not determine (G,P) uniquely, but does determine 
a set ^ such that (G,P) G ^ and such that any e ^ 

could a priori be (G,P). This mathematical framework, called 
optimal uncertainty quantification (OUQ), is based on the 
observation that, given a set of assumptions and information 



about the problem, there exist optimal bounds on uncertainties; 
these are obtained as extreme values of well-defined optimiza- 
tion problems corresponding to extremizing probabilities of 
failure, or of deviations, over the feasible set s^. Observe 
that this framework does not implicitly impose inappropriate 
assumptions, nor does it repudiate relevant information. In- 
deed, as demonstrated in (10 and 11) for the CoM approach, 
OUQ can pose a problem that incorporates the assumptions 
utilized in other common UQ methods (such as Bayesian 
inference [LJH99]) and provide a rigorous optimal bound on 
the uncertainties. 

Although some OUQ problems can be solved analytically, 
most must be solved numerically. To that end, the reduction 
theorems of [OSS 11] reduce the infinite-dimensional feasible 
set to a finite-dimensional subset s^a that has the key 
property that the objective function (PoF) has the same lower 
and upper extreme values over ^ as over 

For example, the reduction for ^md in (10) is to pass to 
measures /i = /ii (8) • • • (8) /i^v such that each marginal measure 
{li is supported on at most two points of the parameter space 
Ei, i.e. {li is a convex combination of two Dirac measures 
(point masses). Having reduced the set of feasible measures 
/i, the set of feasible response functions / is also reduced, 
since we only care about the values of / on the finite support 
of fi and nowhere else. 

We refer the reader to [OSS 11] for the more general 
reduction theorems. The essential point is that if the infor- 
mation/constraints take the form of rii inequalities of the form 
E^- < (for some test functions 0j) and n' inequalities of 
the form [(j)j] < 0, then it is enough to consider {if with 
support on 1 -\-ni-\-n^ points of Ei. 

The reduction theorems leave us with a finite-dimensional 
optimization problem in which the optimization variables are 
suitable parametrizations of the reduced feasible scenarios 
(/,M)- 

A Highly-Configurable Optimization Framework 

We have built a robust optimization framework (mystic) 
[MHA09] that incorporates the mathematical framework de- 
scribed in [OSS 11], and have provided an interface to pre- 
diction, certification, and validation as a framework service. 
The mystic framework provides a collection of optimization 
algorithms and tools that lowers the barrier to solving com- 
plex optimization problems, mystic provides a selection of 
optimizers, both global and local, including several gradient 
solvers. A unique and powerful feature of the framework is 
the ability to apply and configure solver-independent termi- 
nation conditions — a capability that greatly increases the 
flexibility for numerically solving problems with non-standard 
convergence profiles. All of mystic's solvers conform to a 
solver API, thus also have common method calls to configure 
and launch an optimization job. This allows any of mystic's 
solvers to be easily swapped without the user having to write 
any new code. 

The minimal solver interface: 

# the function to be minimized and the initial values 
from mystic .models import rosen as my_moclel 
xO = [0.8, 1.2, 0.7] 
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Figure 1: Conceptual diagram for an optimizer The cost function 
provides a difference metric that accepts input parameters x and 
produces a cost E. 



# configure the solver and obtain the solution 
from mystic . solvers import f min 

solution = fmin (my_model, xO) 

The criteria for when and how an optimization terminates are 
of paramount importance in traversing a function's potential 
well. Standard optimization packages provide a single conver- 
gence condition for each optimizer, mystic provides a set of 
fully customizable termination conditions, allowing the user to 
discover how to better navigate the optimizer through difficult 
terrain. Optimizers can be further configured through several 
methods for choosing the InitialPoints. 
The expanded solver interface: 

# the function to he minimized and initial values 
from mystic .models import rosen as my_model 

xO = [0.8, 1.2, 0.7] 

# get monitor and termination condition objects 
from mystic .monitors import Monitor, VerboseMonitor 
stepmon = VerboseMonitor ( 5 ) 

evalmon = Monitor () 

from mystic .termination import ChangeOverGeneration 
COG = ChangeOverGeneration ( ) 

# instantiate and configure the solver 

from mystic . solvers import NelderMeadSimplexSolver 

solver = NelderMeadSimplexSolver ( len (xO ) ) 

solver . Set InitialPoints (xO) 

solver , SetGenerationMonitor ( stepmon) 

solver . SetEvaluationMonitor (evalmon) 

solver , Solve (my _model, COG) 

# obtain the solution 
solution = solver .bestSolution 

# obtain diagnostic information 

f unction_evals = solver . evaluations 
iterations = solver . generations 
cost = solver , bestEnergy 

# modify the solver configuration, and continue 
COG = ChangeOverGeneration (tolerance^le-8 ) 
solver , Solve (my _model, COG) 

# obtain the new solution 
solution = solver .bestSolution 

mystic provides progress monitors that can be attached 
to an optimizer to track progress of the fitted parameters 
and the value of the cost function. Additionally, monitors 
can be customized to track the function gradient or other 
progress metrics. Monitors can also be configured to record 
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Figure 2: Optimization analysis viewers available in mystic. 
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Figure 3: Basic components provided in the optimizer toolkit. Several 
wrapper classes are also provided for binding components, while 
factory classes are provided for generating components. 



either function evaluations or optimization iterations (i.e. 
generations). For example, using VerboseMonitor ( 5 ) 
in the SetGenerationMonitor method will print the 
bestEnergy to stdout every five generations. 

Constraints Toolkit 

mystic provides a method to constrain optimization to be 
within an A/^-dimensional box on input space, and also a 
method to impose user-defined parameter constraint functions 
on any cost function. Thus, both bounds constraints and 
parameter constraints can be generically applied to any of 
mystic's unconstrained optimization algorithms. Tradition- 
ally, constrained optimization problems tend to be solved 
iteratively, where a penalty is applied to candidate solutions 
that violate the constraints. Decoupling the solving of con- 
straints from the optimization problem can greatly increase the 
efficiency in solving highly-constrained nonlinear problems ~ 
effectively, the optimization algorithm only selects points that 
satisfy the constraints. Constraints can be solved numerically 
or algebraically, where the solving of constraints can itself 
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be cast as an optimization problem. Constraints can also be 
dynamically applied, thus altering an optimization in progress. 

Penalty-based methods indirectly modify the candidate so- 
lution by applying a change in energy AE" = k- p{x) to the 
unconstrained cost function f{x) when the constraints are 
violated. The modified cost function (j) is thus written as: 



(l){x)=f{x)^k-p{x) 



(13) 



Set-based methods directly modify the candidate solution by 
applying a constraints solver c that ensures the optimizer 
will always select from a set of candidates that satisfy the 
constraints. The constraints solver has an interface x' = c{x), 
and the cost function becomes: 



0(x)=/(c(x)) 



(14) 



Adding parameter constraints to a solver is as simple as build- 
ing a constraints function, and using the SetConstraint s 
method. Additionally, simple bounds constraints can also be 
applied through the Set Strict Ranges method: 

# a user-provided constraints function 
def constrain (x) : 

x[l] = x[0] 
return x 

# the function to be minimized and the bounds 
from mystic .models import rosen as my_model 
lb = [0.0, 0.0, 0.0] 

ub = [2.0, 2.0, 2.0] 

# get termination condition object 

from mystic .termination import ClnangeOverGeneration 
COG = CliangeOverGeneration ( ) 

# instantiate and configure the solver 

from mystic . solvers import NelderMeadSimplexSolver 

solver = NelderMeadSimplexSolver ( len (xO ) ) 

solver . SetRandomlnitialPoints (lb, ub) 

solver . SetStrictRanges (lb, ub) 

solver . SetConstraint s (constrain) 

solver . Solve (my _model, COG) 

# obtain the solution 
solution = solver . bestSolution 

mystic provides a simple interface to a lot of underlying 
complexity ~ thus allowing a non-specialist user to easily 
access optimizer configurability and high-performance com- 
puting without a steep learning curve. This feature must also 
be applied to the application of constraints on a function or 
measure. The natural syntax for a constraint is one of symbolic 
math, hence mystic leverages SymPy [SYMll] to construct 
a symbolic math parser for the translation of the user's input 
into functioning constraint code objects: 

# a user-provided constraints function 

constraints = """ 

x2 = xl 
II II II 

from mystic . constraints import parse 
constrain = parse ( constraint s ) 

The constraints parser is a constraints factory method that can 
parse multiple and nonlinear constraints, hard or soft (i.e. "~") 
constraints, and equality or inequality (i.e. ">") constraints. 

Similar tools exist for creating penalty functions, including a 
SetPenalty method for solvers. Available penalty methods 
include the exterior penalty function method [VEN09], the 
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Figure 4: Conceptual diagram for a service-based model. Here, the 
job is the fundamental commodity of work, and is the object on which 
the service is based - in mystic, this is typically the user's model or 
a cost function. Services have a global unique identifier, and thus can 
easily be called by proxy. Note that services may not be located on 
the machine that requested the service be spawned. Services also can 
be imbued with infrastructure for monitoring and handling events. 
Monitors write to a stream that can be piped into another object, 
such as a logger or one of mystic's viewers. 



augmented Lagrange multiplier method [KSK94], and the 
logarithmic barrier method [JJB03]. At the low-level, penalty 
functions are bound to the cost function using mystic's 
f unct ionWrapper method. 

It is worth noting that the use of a constraints solver 
c does not require the constraints be bound to the cost 
function. The evaluation of the constraints are decoupled from 
the evaluation of the cost function ~ hence, with mystic, 
highly-constrained optimization decomposes to the solving 
of K independent constraints, followed by an unconstrained 
optimization over only the set of valid points. This method 
has been shown effective for solving optimization problems 
where ^^200 [OSS 11]. 

Seamless Migration to Parallel Computing 

mystic is built from the ground up to utilize parallel and 
distributed computing. The decomposition of optimization 
algorithms into their component parts allow this decomposition 
to not only be in an abstraction layer, but across process- 
space, mystic provides a mode IFact cry method that 
convers a user's model to a service. We define a service to 
be an entity that is callable by globally unique identifier. 
Services can also be called by proxy. In mystic, services 
also include infrastructure for monitoring and handling events. 
An optimization is then composed as a network of interacting 
services, with the most common being the user's model or 
cost function being mapped over parallel resources. 

mystic provides several stock models and model factories 
that are useful for testing: 

# generate a model from a stock ^model factory'' 
from mystic .models . lorentzian import Lorentzian 
lorentz = Lorent zian ( coef f s ) 

# evaluate the model 
y = lorentz (x) 

Model factory methods insert pathos infrastructure, thus 
casting a model as a callable service that has been imbued with 
pathos infrastructure as shown in Figure (4). The default 
launcher and map included in mystic are functionally 
equivalent to execution and map within the standard Python 
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Figure 5: Use of a modelFactory to cast a user's model F(x) as 
a service. The model and experimental data G are then bound with a 
costFactory to produce a cost function. A costFactory can 
accept a raw user's model, a model proxy, or a model service (as 
shown here). A typical metric is \F(x) — G\^. 

distribution. Any user-provided function can be cast as a 
service through the use of a modelFactory: 

# a user-provided model function 
def identify (x) 

return x 

# add pathos infrastructure (included in mystic) 
from mystic . tools import modelFactory, Monitor 
evalmon = Monitor () 

my_model = modelFactory ( identify , monitor=evalmon) 

# evaluate the model 
y = my_model (x) 

# evaluate the model with a map function 
from mystic. tools import PythonMap 
my_map = PythonMap ( ) 

z = my_map (my_model, range (10)) 

A Framework for Heterogeneous Computing 

We have developed a computational job management frame- 
work (pathos) [MBA 10] that offers a simple, efficient, 
and consistent user experience in a variety of heterogeneous 
environments from multi-core workstations to networks of 
large-scale computer clusters, pathos provides a single en- 
vironment for developing and testing algorithms locally ~ and 
enables the user to execute the algorithms on remote clusters, 
while providing the user with full access to their job history, 
pathos primarily provides the communication mechanisms 
for configuring and launching parallel computations across 
heterogenous resources, pathos provides stagers and launch- 
ers for parallel and distributed computing, where each launcher 
contains the syntactic logic to configure and launch jobs in an 
execution environment. Some examples of included launchers 
are: a queue-less MPI-based launcher, a SSH-based launcher, 
and a multiprocessing launcher, pathos also provides 
a map-reduce algorithm for each of the available launchers, 
thus greatly lowering the barrier for users to extend their code 
to parallel and distributed resources, pathos provides the 
ability to interact with batch schedulers and queuing systems, 
thus allowing large computations to be easily launched on 
high-performance computing resources. One of the most pow- 
erful features of pathos is sshTunnel, which enables a 
user to automatically wrap any distributed service calls within 
an SSH tunnel. 

pathos is divided into four subpackages: dill (a util- 
ity for serialization of Python objects), pox (utilities for 
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filesystem exploration and automated builds), pyina (a MPI- 
based parallel mapper and launcher), and pathos (distributed 
parallel map-reduce and SSH communication). 

pathos utilizes pyre, which provides tools for connecting 
components and managing their interactions. The core com- 
ponent used by pathos is a service ~ a callable object with 
a configurable connection mechanism. A service can utilize 
Launcher and Monitor objects (which provide abstraction 
to execution and logging, respectively), as well as Strategy 
objects (which provide abstraction patterns for coupling ser- 
vices). A standard interface for services enables massively 
parallel applications that utilize distributed resources to be 
constructed from a few simple building blocks. A Launcher 
contains the logic required to initiate execution on the cur- 
rent execution environment. The selection of launcher will 
determine if the code is submitted to a batch queue, run 
across SSH tunneled RPC connections, or run with MPI on 
a multiprocessor. A Strategy provides an algorithm to 
distribute the workload among available resources. Strategies 
can be static or dynamic. Examples of static strategies include 
theequalportion strategy and the carddealer strategy. 
Dynamic strategies are based on the concept of a worker 
pool, where there are several workload balancing options to 
choose from. Strategies and launchers can be coupled together 
to provide higher-level batch and parallel-map algorithms. A 
Map interface allows batch processing to be decoupled from 
code execution details on the selected platforms, thus enabling 
the same application to be utilized for sequential, parallel, and 
distributed parallel calculations. 

Globally Unique Message Passing 

We must design for the case where an optimizer's calculation 
spans multiple clusters, with a longevity that may exceed the 
uptime of any single cluster or node, hydra enables any 
Python object to obtain a network address. After obtaining 
an address, an object can asynchronously exchange messages 
with other objects on the network. Through the use of proxy 
objects, sending messages to remote objects is easy as calling 
an instance method on a local object. A call to a proxy trans- 
parently pickles the function name along with the arguments, 
packages the message as a datagram, and sends it over the 
network to the remote object represented by the proxy. On 
the recieving end, there is a mechanism for responding to 
the sender of the current message. Since message sending is 
asynchronous, an object responds to a message by sending 
another message. 

The modelFactory method essentially provides mystic 
with a high-level interface for a pathos server, with an 
option to bind a monitor directly to the service. The lower- 
level construction of a distributed service, using SSH-based 
communication, is as follows: 

# a user-provided model function 
def identify (x) 

return x 

# cast the model as a distributed service 
from pathos . servers import sshServer 

id = 'foo.caltech.edu:50000:spike42' 
my_proxy = sshServer ( identify , server=id) 
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Figure 6: Conceptual diagram for heterogeneous job management. 
A distributed parallel map function is used to copy a service n times 
on N machines. If the object being mapped is not a service, then 
the services manager is omitted from the diagram - the jobs still 
undergo a distributed launch, but are managed at the machine level. 



# evaluate the model via proxy 
y = my_proxy(x) 

Parallel map functions are built around available launchers, 
providing a high-level interface to launching several copies of 
a model in parallel. The creation of a parallel map that will 
draw from a pool of two local workers and all available IPC 
servers at ' f oo . caltech . edu' is shown below: 

# a user-provided model function 
def identify (x) 

return x 

# select and configure a parallel map 
from pathos. maps import ipcPool 

my_map = ipcPool(2, servers= foo . caltech . eduM ) 

# evaluate the model in parallel 
z = my_map ( identify , range (10)) 



Serialization 

dill extends Python's pickle module for serializing and 
de- serializing Python objects to the majority of the built- 
in Python and NumPy types. Serialization is the process of 
converting an object to a byte stream, the inverse of which is 
converting a byte stream back to a Python object hierarchy. 

dill provides the user the same interface as the pickle 
module, and also includes some additional features. In addition 
to pickling Python objects, dill provides the ability to save 
the state of an interpreter session in a single command. Hence, 
it would be feasible to save a interpreter session, close the 
interpreter, ship the pickled file to another computer, open a 
new interpreter, unpickle the session and thus continue from 
the "saved" state of the original interpreter session. 

Filesystem Interaction 

pox provides a collection of utilities for navigating and 
manipulating filesy stems. This module is designed to facilitate 



some of the low level operating system interactions that are 
useful when exploring a filesystem on a remote host, where 
queries such as "what is the root of the filesystem?", "what is 
the user's name?", and "what login shell is preferred?" become 
essential in allowing a remote user to function as if they were 
logged in locally. While pox is in the same vein of both the os 
and shutil built-in modules, the majority of its functionality 
is unique and compliments these two modules. 

pox provides Python equivalents of several unix shell 
commands such as "which" and "find". These commands 
allow automated discovery of what has been installed on an 
operating system, and where the essential tools are located. 
This capability is useful not only for exploring remote hosts, 
but also locally as a helper utility for automated build and 
installation. 

Several high-level operations on files and filesystems are 
also provided. Examples of which are: finding the location 
of an installed Python package, determining if and where the 
source code resides on the filesystem, and determining what 
version the installed package is. 

pox also provides utilities to enable the abstraction of 
commands sent to a remote filesystem. In conjunction with 
a registry of environment variables and installed utilites, pox 
enables the user to interact with a remote filesystem as if they 
were logged in locally. 

Distributed Staging and Launching 

pathos provides methods for configuring, launching, moni- 
toring, and controlling a service on a remote host. One of the 
most basic features of pathos is the ability to configure and 
launch a IPC-based service on a remote host, pathos seeds 
the remote host with a small portpicker script, which 
allows the remote host to inform the localhost of a port that 
is available for communication. 

Beyond the ability to establish a IPC service, and then 
post requests, is the ability to launch code in parallel. Unlike 
parallel computing performed at the node level (typically with 
MPI), pathos enables the user to launch jobs in parallel 
across heterogeneous distributed resources, pathos provides 
a distributed map-reduce algorithm, where a mix of local pro- 
cessors and distributed IPC services can be selected, pathos 
also provides a very basic automated load balancing service, as 
well as the ability for the user to directly select the resources. 

A high-level interface is provided which yields a map- 
reduce implementation that hides the IPC internals from the 
user. For example, with ipcPool, the user can launch their 
code as a distributed parallel service, using standard Python 
and without writing a line of server or parallel batch code, 
pathos also provides tools to build a custom Map. In 
following code, the map is configured to 'autodetect' 
the number of processors, and only run on the localhost: 

# configure and build map 

from pathos . launchers import ipc 

from pathos . strategies import pool 

from pathos. tools import mapFactory 

my_map = mapFactory ( launcher=ipc, st rategy=pool ) 

IPC servers and communication in general is known to be 
insecure. However, instead of attempting to make the IPC 
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communication itself secure, pathos provides the ability to 
automatically wrap any distributes service or communication 
in an SSH tunnel. SSH is a universally trusted method. 
Using sshTunnel, pathos has launched several distributed 
calculations on clusters at National Laboratories, and to date 
has performed test calculations that utilize node-to-node com- 
munication between two national lab clusters and a user's 
laptop, pathos allows the user to configure and launch 
at a very atomistic level, through raw access to ssh and 
sop. Any distributed service can be tunneled, therefore less- 
secure methods of communication can be provided with secure 
authentication: 

# establish a tunnel 

from pathos . tunnel import sshTunnel 
uid = ' foo . caltech . edu : 12345 : tunnel69' 
tunnel_proxy = sshTunnel (uid) 

# inspect the ports 
localport = tunnel_proxy . Iport 
remoteport = tunnel_proxy . rport 

# a user-provided model function 
def identify (x) 

return x 

# cast the model as a distributed service 
from pathos . servers import ipcServer 

id = ' localhost : %.s:bug01' % localport 
my_proxy = ipcServer ( identify , server=id) 

# evaluate the model via tunneled proxy 
y = my_proxy(x) 

# disconnect the tunnel 
tunnel_proxy . disconnect ( ) 

Parallel Staging and Launching 

The pyina package provides several basic tools to make MPI- 
based high-performance computing more accessable to the end 
user. The goal of pyina is to allow the user to extend their 
own code to MPI-based high-performance computing with 
minimal refactoring. 

The central element of pyina is the parallel map-reduce 
algorithm, pyina currently provides two strategies for exe- 
cuting the parallel-map, where a strategy is the algorithm for 
distributing the work list of jobs across the availble nodes. 
These strategies can be used "in-the-raw" (i.e. directly) to 
provide map-reduce to a user's own MPI-aware code. Further, 
pyina provides several map-reduce implementations that 
hide the MPI internals from the user. With these Map objects, 
the user can launch their code in parallel batch mode ~ using 
standard Python and without ever having to write a line of 
Parallel Python or MPI code. 

There are several ways that a user would typically launch 
their code in parallel ~ directly with mpirun ormpiexec, 
or through the use of a scheduler such as torque or slurm. 
pyina encapsulates several of these launching mechanisms 
as Launchers, and provides a common interface to the 
different methods of launching a MPI job. In the following 
code, a custom Map is built to execute MPI locally (i.e. not 
to a scheduler) using the carddealer strategy: 

# configure and build map 

from pyina . launchers import mpirun 
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Figure 7: Conceptual diagram for a carddealer-DE optimizer 
The optimizer contains a map function that stages n copies of the 
user's model F{x) in parallel across distributed resources. 



from pyina . strategies import carddealer as card 
from pyina. tools import mapFactory 

my_map = mapFactory ( 4 , launcher=mpirun, st rategy=card) 

New Massively-Parallel Optimization Algorithms 

In mystic, optimizers have been extended to parallel when- 
ever possible. To have an optimizer execute in parallel, 
the user only needs to provide the solver with a paral- 
lel map. For example, extending the Differential Evolution 
[SKP95] solver to parallel is involves passing a Map to the 
SetEvaluationMap method. In the example below, each 
generation has 20 candidates, and will execute in parallel using 
MPI with 4 workers: 

# the function to be minimized and the bounds 
from mystic .models import rosen as my_model 
lb = [0.0, 0.0, 0.0] 

ub = [2.0, 2.0, 2.0] 

# get termination condition object 

from mystic .termination import CliangeOverGeneration 
COG = CliangeOverGeneration ( ) 

# select the parallel launch configuration 
from pyina. maps import MpirunCarddealer 
my_map = MpirunCarddealer ( 4 ) 

# instantiate and configure the solver 

from mystic . solvers import Dif f erentialEvolutionSolver 
solver = Dif f erentialEvolutionSolver ( len ( lb) , 20) 
solver , SetRandomlnitialPoints (lb, ub) 
solver , SetStrictRanges ( lb, ub) 
solver . SetEvaluationMap (my_map) 
solver , Solve (my_model, COG) 

# obtain the solution 
solution = solver . bestSolution 

Another type of new parallel solver utilizes the 
SetNestedSolver method to stage a parallel launch 
of N optimizers, each with different initial conditions. The 
following code shows the BuckshotSolver scheduling a 
launch of A/^ = 20 optimizers in parallel to the default queue, 
where 5 nodes each with 4 processors have been requested: 

# the function to be minimized and the bounds 
from mystic .models import rosen as my_model 
lb = [0.0, 0.0, 0.0] 

ub = [2.0, 2.0, 2.0] 
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Figure 8: Conceptual diagram for a lattice-Powell optimizer 
N Powell's local-search optimizers are launched in parallel, with 
each optimizer starting from the center of a different lattice cuboid 
in parameter space. A buckshot -Powell optimizer is similar; 
however, instead utilizes a uniform random distribution of initial 
values. 



# get monitor and termination condition objects 
from mystic .monitors import LoggingMonitor 
stepmon = LoggingMonitor ( 1 , ^log.txt') 

from mystic . termination import ChangeOverGeneration 
COG = ChangeOverGeneration ( ) 

# select the parallel launch configuration 
from pyina.maps import TorqueMpirunCarddealer 
my_map = TorqueMpirunCarddealer ( ^ 5 : ppn^4 ^ ) 

# instantiate and configure the nested solver 
from mystic . solvers import PowellDirectionalSolver 
my_solver = PowellDirectionalSolver ( len ( lb) ) 
my_solver . SetStrictRanges (lb, ub) 

my_solver . SetEvaluationLimit s ( 50 ) 

# instantiate and configure the outer solver 
from mystic . solvers import BuckshotSolver 
solver = BuckshotSolver ( len ( lb) , 20) 
solver , SetRandomlnitialPoints (lb, ub) 
solver . SetGenerationMonitor ( stepmon) 
solver . SetNestedSolver (my_solver ) 

solver , SetSolverMap (my_map) 
solver . Solve (my _model, COG) 

# obtain the solution 
solution = solver . bestSolution 



Probability and Uncertainty Tooklit 

The software framework presented in this paper was de- 
signed to solve UQ problems. Calculation of the upper and 
lower bounds for probability of failure is provided as a 
framework service. The McDiarmid subdiameter is a model- 
based measure of sensitivity, and is cast within mystic as 
a global optimization. Diameter calculations can be coupled 
with partitioning algorithms, and used to discover regions of 
critical behavior. Optimization over probability measures is 
also available as a framework service, and is utilized in (OUQ) 
calculations of optimal bounds. 

The minimization or maximization of a cost function is 
the basis for performing most calculations in mystic. The 
optimizer generates new trial parameters, which are evaluated 
in a user-provided model function against a user-provided 




diameters help zero-in 
on regions of parameter 
space where parameters 
have the desired impact 



Figure 9: Coupling an iterative partitioning algorithm with a 
sensitivity calculation enables the discovery of critical regions in 
parameter space. 



metric. Two simple difference metrics provided are: metric = 
— Gp, where F is the model function evaluated at some 
trial set of fit parameters and G is the corresponding 
experimental data ~ and metric = \F {x) — F {y)\'^ , where x and 
y are two slightly different sets of input parameters (6). 

mystic provides factory methods to automate the gener- 
ation of a cost function from a user's model. Conceptually, a 
CO St Fact cry is as follows: 

# prepare a (F (X) - G) **2 a metric 
def costFactory (my_model, my_data) : 

def cost (param) : 

# compute the cost 

return ( my_model (param) - my_data ) ^^2 
return cost 

Suboscillations (6), used in calculations of rigorous sensi- 
tivity (such as Di/D), can also be cast as a cost metric: 

# prepare a (F (X) - F (X' ) ) cost metric 
def suboscillationFactory (my_model, i) : 

def cost (param) : 

# get X and X' (Xi' is appended to X at param [-1]) 
X = param [ : -1 ] 

x_prime = param[:i] + param[-l:] + param [ i+1 : -1 ] 

# compute the suboscillation 

return -( my_model (x) - my_model (x_prime) ) ^^2 
return cost 

The diameter D (7) is the root-mean square of its component- 
wise suboscillations. The calculation of the diameter is per- 
formed as a nested optimization, as shown above for the 
BuckshotSolver. Each inner optimization is a calculation 
of a component suboscillation, using the a global optimizer 
(such as Dif f erentialEvolutionSolver) and the cost 
metric shown above. 

The optimization algorithm takes a set of model parameters 
^ and the current measure of oscillation 0{^) as inputs, and 
produces an updated The optimization loop iterates until 
the termination conditions are satisfied. 

When the global optimization terminates the condition 
0{^) < —{oscj-\-£) is satisfied, and the final set ^ is 
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which reduces to the finite-dimensional subset 
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Figure 10: Optimal uncertainty quantification is an optimization 
of probability measures over design parameter space. Collapse of 
probability masses corresponds to the determination of the critical 
design parameters. 



composed of X and X\ 

OUQ problems can be thought of optimization problems 
where the goal is to find the global maximum of a probability 
function ii[H < 0], where < is a failure criterion for 
the model response function H. Additional conditions in an 
OUQ problem are provided as constraints on the information 
set. Typically, a condition such as a mean constraint on H, 
mi < IE^[//] < m2, will be imposed on the maximization. 
After casting the OUQ problem in terms of optimization and 
constraints, we can plug these terms into the infrastructure 
provided by my s t i c . 

Optimal uncertainty quantification (OUQ) is maximization 
over a probability distribution, and not over a standard dif- 
ference metric. Therefore the fundamental data structure is 
not the user-provided model function, but is a user-configured 
probability measure. For example, a discrete measure is rep- 
resented by a collection of support points, each with an 
accompanying weight. Measures come with built-in methods 
for calculating the mass, range, and mean of the measure, and 
also for imposing a mass, range, and mean on the measure. 
Measures also have some very basic operations, including 
point addition and subtraction, and the formation of product 
measures. 

Global optimizations used in solving OUQ problems are 
composed in the same manner as shown above for the 
Dif f erentialEvolutionSolver. The cost function, 
however, is not formulated as in the examples above ~ OUQ 
is an optimization over product measures, and thus uses 
mystic's product_measure class as the target of the 
optimization. Also as shown above, the bounds constraints 
are imposed with the Set Strict Ranges method, while 
parameter constraints (composed as below) are imposed with 
the SetConstraints method. The union set of these con- 
straints defines the set 

So for example, let us define the feasable set 



/ = iny_inodel : flLi [lt>i^^t)i] 
Al = (8)-=i e (8)-=i ^([lbi,ubi]), 



for X and y e Y[^=i[^'b±,uh±\,' 
and w e [0,1], 



(16) 



(15) 



where x = some (xi,X2,X3), y = some (ji, 3^27^3), and w = 
some (wi,W2,W3). 

To solve this OUQ problem, we first write the code for the 
bounds, cost function, and constraints ~ then we plug this 
code into a global optimization script, as noted above. 

OUQ requires the user provide a list of bounds 
that follow the formatting convention that mystic's 
product_measure .load uses to build a product measure 
from a list of input parameters. This roughly follows the 
definition of a product measure as shown in equation (16), 
and also is detailed in the comment block below: 

# OUQ requires bounds in a very specific form. . . 

# param = [wxi ] *nx + [xi ] *nx \ 

# + [wyi ] y^rny + [yi ] *ny \ 

# + [wzi ] *nz + [zi ] *nz 
npts = (nx,ny,nz) 

lb = (nx ^ w_lower) + (nx ^ x_lower) \ 

+ (ny * w_lower) + (ny ^ y_lower) \ 

+ (nz ^ w_lower) + (nz ^ z_lower) 

ub = (nx * w_upper) + (nx ^ x_upper) \ 

+ (ny * w_upper) + (ny ^ y_upper) \ 

+ (nz ^ w_upper) + (nz ^ z_upper) 

The constraints function and the cost function typically require 
the use of measure mathematics. In the example below, the 
constraints check if measure. mass ~ 1.0; if not, the the 
measure's mass is normalized to 1.0. The second block of 
constraints below check if mi < IE^[//] < m2, where mi = 
target_mean — error and m2 = target_mean + 
error; if not, an optimization is performed to satisfy this 
mean constraint. The product_measure is built (with 
load) from the optimization parameters param, and after 
all the constraints are applied, flatten is used to extract 
the updated param: 

from mystic .math .measures import split_param 

from mystic .math . dirac_measure import product_measure 

from mystic. math import almostEqual 

# split bounds into weight-only & sample-only 
w_lb, m_lb = split_param ( lb, npts) 

w_ub, m_ub = split_param (ub, npts) 

# generate constraints function 
def constraints (param) : 

prodmeasure = product_measure ( ) 
prodmeasure . load (param, npts) 

# impose norm on measures 
for measure in prodmeasure: 

if not almostEqual ( float (measure . mass ) , 1.0): 
measure . normalize ( ) 

# impose expectation on product measure 

E = float (prodmeasure . get_expect (my_model) ) 
if not (E <= float (target_mean + error) ) \ 
or not ( float (target_mean - error) <= E) : 

prodmeasure , set_expect ( (target_mean, error), \ 

my_model, (m_lb, m_ub) ) 

# extract weights and positions 
return prodmeasure . flatten ( ) 
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The PoF is calculated in the cost function with the pof 
method: 

# generate maximizing function 

def cost (param) : 

prodmeasure = product_measure ( ) 

prodmeasure . load (param, npts) 

return MINMAX ^ prodmeasure . pof (my_model ) 

We find the supremum (as in 12) when MINMAX=-1 and, upon 
solution, the function maximum is -solver.bestEnergy. 
We find the infimum when MINMAX=1 and, upon solution, the 
function minimum is solver.bestEnergy. 

Future Developments 

Many of the features presented above are not currently in 
released versions of the code. Of primary importance is to 
migrate these features from development branches to a new 
release. 

The next natural question beyond "what is the sensitivity of 
a model to an input parameter?" is "how does the correlation 
between input parameters affect the outcome of the model?". 
Methods for calculating parameter correlation will be very 
useful in analysis of results. Another natural question is how 
to handle uncertainty in the data. 

New partitioning algorithms for the discovery of regions of 
critical behavior will be added to mystic. Currently the only 
partitioning rule drives the optimizer toward partitioning space 
such that the upper bounds of a "piecewise-McDiarmid" type 
are iteratively tightened [STMll]. We will extend the parti- 
tioning algorithm not to refine the diameter, but to discover 
regions where the diameters meet a set of criteria (such as: 
regions where there are two roughly equal subdiameters that 
account for 90% or more of the total diameter (i.e. automated 
discovery of regions where two parameters compete to govern 
the system behavior), mystic will also further expand its 
base of available statistical and measure methods, equation 
solvers, and also make available several more traditional 
uncertainty quantification methods, mystic will continue to 
expand its base of optimizers, with particular emphasis on 
new optimization algorithms that efficiently utilize parallel 
computing, mystic currently has a few simple parallel 
optimization algorithms, such as the LatticeSolver and 
BuckshotSolver solvers; however, algorithms that utilize a 
variant of game theory to do speculation about future iterations 
(i.e. break the paradigm of an iteration being a blocker to 
parallelism), or use parallelism and dynamic constraints to 
allow optimizers launched in parallel to avoid finding the 
same minimum twice, are planned. Parallelism in optimization 
also allows us to not only find the global minima, but to 
simultaneously find all local minima and transition points - 
- thus providing a much more efficient means of mapping out 
a potential energy surface. Solving uncertainty quantification 
problems requires a lot of computational resources and often 
must require a minimum of both model evaluations and 
accompanying experiments, so we also have to keep an eye 
on developing parallel algorithms for global optimization with 
overall computational efficiency. 

pathos includes utilities for filesystem exploration and 
automated builds, and a utility for the serialization of Python 



objects, however these framework services will need to be 
made more robust as more platforms and more extensive 
objects and codes are tackled. Effort will continue on expand- 
ing the management and platform capabilities for pathos, 
unifying and hardening the map interface and providing load 
balancing for all types of connections. The high-level inter- 
face to analysis circuits will be extended to encompass new 
types of logic for combining and nesting components (as 
nested optimizers are utilized in many materials theory codes). 
Monitoring and logging to files and databases across parallel 
and distributed resources will be migrated from mystic and 
added as pathos framework services. 

Summary 

A brief overview of the mathematical and software compo- 
nents used in building a software framework for predictive 
science is presented. 
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